#雷电波形分解，重点分析峰值、上升时间、电流/电压积分和电流电压平方的积分，以及频率分量
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simpson
from scipy.signal import butter, filtfilt, freqz
import pandas as pd


#数据前处理，将波形起点置为0
#原始数据，最初测试波形，实际测试波形峰值，在实际的过程中一般应该是1：1的
def data_prepocessing(data,init_level,actual_eval):
    #对于波形数据必须从(0,0)开始,因此先将数据小于0的部分去掉
    time=data[:,0]*1e6
    ratio=actual_eval/init_level
    amplitude=data[:,1]
    # plt.figure(1)
    # plt.rcParams['font.sans-serif'] = ['SimHei'] 
    # plt.plot(time,amplitude)
    # plt.xlabel("时间/us")
    # plt.ylabel("幅值/V")
    # plt.title("测试所得原始波形")
    # 去掉时间小于0的部分
    index = 0
    n = len(time)
    while index < n and time[index] <= 0:
        index += 1
    # 然后截取数组
    time = time[index:]
    print(f"当前信号起点{time[0]}")
    amplitude = amplitude[index:]*ratio
    print(f"当前信号起点处的幅值{amplitude[0]}")

    #此时得到的值
    dt=time[0]
    print(f"当前信号起点{dt}")
    # da=amplitude[0]
    time=time-dt
    # amplitude=amplitude-da
    # plt.figure(2)
    # plt.rcParams['font.sans-serif'] = ['SimHei'] 
    # plt.plot(time,amplitude)
    # plt.xlabel("时间/us")
    # plt.ylabel("幅值/V")
    # plt.title("雷电信号起点归0后的波形")
    peak_zero=max(amplitude)
    print(f"归零后的波形峰值为{peak_zero}V")
    plt.rcParams['axes.unicode_minus']=False 
    plt.show()
    return time,amplitude

#提取当前波形的参数
def lightning_abstract_separate(time,amplitude):#假设当前的雷电时域数据为us和A
    time=time/(1e6)
    
    amplitude=amplitude
    # plt.figure(3)
    # plt.plot(time,amplitude)
    # plt.xlabel("时间/s")
    # plt.ylabel("幅值/V")
    # plt.title("雷电原始波形")
    
    #信号峰值
    peak=max(amplitude)
    print(f"波形峰值为{peak}V")


    #提取上升沿的上升时间10%到90%的时间
    for i in range(len(time)):
        if (abs(amplitude[i]-peak*0.1)/(peak*0.1))<0.01:
            peak_forward_time=float(time[i])
            peak_forward_amplitude=float(amplitude[i])
            print(f"当前信号的10%峰值时间为{peak_forward_time},当前峰值为{peak_forward_amplitude}")
            break
    for i in range(len(time)):
        if (abs(amplitude[i]-peak*0.9)/(peak*0.9))<0.01:
            peak_backforward_time=float(time[i])
            peak_backforward_amplitude=float(amplitude[i])
            print(f"当前信号的90%峰值时间为{peak_backforward_time},当前峰值为{peak_backforward_amplitude}")
            break
    rate_of_rise=(peak_backforward_amplitude-peak_forward_amplitude)/((peak_backforward_time-peak_forward_time))
    print(f"波形上升速率为{rate_of_rise}A/s")


    #计算电荷
    content=np.trapz(amplitude,time)
    print(f"当前的电荷量为{content}V-s")

    #计算能量
    amplitude_power=pow(amplitude,2)
    energy=np.trapz(amplitude_power,time)
    print(f"当前的能量为{energy}V^2s")


    #计算数据在1MHz的分量
    fs=10e6
    lowcut=990e3
    highcut=1010e3
    order=6
    nyq = 0.5 * fs  # 奈奎斯特频率
    low = lowcut / nyq
    high = highcut / nyq
    amplitude_1MHz=butter_bandpass_filter(amplitude,lowcut,highcut,fs,order)
    b, a = butter_bandpass(lowcut, highcut, fs, order)
    w, h = freqz(b, a, worN=2000)
    print(f"当前最大的响应为{max(abs(h))}")
    #考虑滤波器的响应
    amplitude_1MHz=amplitude_1MHz/(max(abs(h)))
    # plt.figure(4)
    # plt.plot(time,amplitude_1MHz)
    # plt.xlabel("时间/s")
    # plt.ylabel("幅值/A")
    # plt.title("雷电信号完成滤波之后的波形")
    # plt.show()
    peak_1MHz=max(amplitude_1MHz)
    print(f"1MHz波形峰值为{peak_1MHz}A")
    # plt.show()
    #完成滤波的波形，提取上升时间
    for i in range(len(time)):
        if (abs(amplitude_1MHz[i]-peak_1MHz*0.1)/(peak_1MHz*0.1))<0.01:
            peak_forward_time=float(time[i])
            peak_forward_amplitude=float(amplitude_1MHz[i])
            print(f"当前信号的10%峰值时间为{peak_forward_time},当前峰值为{peak_forward_amplitude}")
            break
    for i in range(len(time)):
        if (abs(amplitude_1MHz[i]-peak_1MHz*0.9)/(peak_1MHz*0.9))<0.01:
            peak_backforward_time=float(time[i])
            peak_backforward_amplitude=float(amplitude_1MHz[i])
            print(f"当前信号的90%峰值时间为{peak_backforward_time},当前峰值为{peak_backforward_amplitude}")
            break
    rate_of_rise_1MHz=(peak_backforward_amplitude-peak_forward_amplitude)/(peak_backforward_time-peak_forward_time)
    print(f"1MHz波形上升速率为{rate_of_rise_1MHz}A/s")
    return peak,content,energy,rate_of_rise,rate_of_rise_1MHz


#定义滤波器
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs  # 奈奎斯特频率
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    # 使用filtfilt实现零相位滤波（避免相位偏移）
    y = filtfilt(b, a, data)
    return y

#确定当前的波形组合
def wave_combind_voltage(peak,content,energy,rate_of_rise,rate_of_rise_1MHz):
    num_list=[["2波1电平",50,1.18e-4,1.96e-3,2.73e8],
              ["2波2电平",125,2.94e-4,1.22e-2,6.83e8],
              ["2波3电平",300,7.07e-4,7.05e-2,1.64e9],
              ["2波4电平",750,1.77e-3,4.41e-1,4.1e9],
              ["2波5电平",1600,3.77e-3,2,8.74e9],

              ["4波1电平",50,4.73e-3,1.25e-1,3.5e7],
              ["4波2电平",125,1.18e-2,7.81e-1,8.75e7],
              ["4波3电平",300,2.84e-2,4.5,2.1e8],
              ["4波4电平",750,7.09e-2,28.1,5.25e8],
              ["4波5电平",1600,1.51e-1,128,1.12e8],

              ["3波1电平",100,0,0,6.65e8],
              ["3波2电平",250,0,0,1.66e9],
              ["3波3电平",600,0,0,3.99e9],
              ["3波4电平",1500,0,0,9.98e9],
              ["3波5电平",3200,0,0,2.13e10],

              
              ]
    
    # num_list_2wave=[
    #           ["2波1电平",50,1.18e-4,1.96e-3,2.73e8],
    #           ["2波2电平",125,2.94e-4,1.22e-2,6.83e8],
    #           ["2波3电平",300,7.07e-4,7.05e-2,1.64e9],
    #           ["2波4电平",750,1.77e-3,4.41e-1,4.1e9],
    #           ["2波5电平",1600,3.77e-3,2,8.74e9],
    # ]
    # num_list_4wave=[
    #           ["4波1电平",50,4.73e-3,1.25e-1,3.5e7],
    #           ["4波2电平",125,1.18e-2,7.81e-1,8.75e7],
    #           ["4波3电平",300,2.84e-2,4.5,2.1e8],
    #           ["4波4电平",750,7.09e-2,28.1,5.25e8],
    #           ["4波5电平",1600,1.51e-1,128,1.12e8],
    # ]
    # num_list_3wave=[
    #           ["3波1电平",100,0,0,6.65e8],
    #           ["3波2电平",250,0,0,1.66e9],
    #           ["3波3电平",600,0,0,3.99e9],
    #           ["3波4电平",1500,0,0,9.98e9],
    #           ["3波5电平",3200,0,0,2.13e10],
    # ]
    evaluate=np.zeros([15,4])
    for i in range(len(num_list)):
        if num_list[i][1]>peak:
            evaluate[i][0]=True

    for i in range(len(num_list)):
        if num_list[i][2]>content:
            evaluate[i][1]=True

    for i in range(len(num_list)):
        if num_list[i][3]>energy:
            evaluate[i][2]=True

    for i in range(len(num_list)-5):
        if num_list[i][4]>rate_of_rise:
            evaluate[i][3]=True

    for i in range(len(num_list)-10):
        if num_list[i+10][4]>rate_of_rise_1MHz:
            evaluate[i+10][3]=True
    print(evaluate)
    wave_single=[["wave_name_1&level"]]
    for row in range(len(evaluate)):
        if np.all(evaluate[row]):
            print(f"当前波形可由{num_list[row][0]}描述")
            wave_single.append({num_list[row][0]})

    print(wave_single)
    file = open(r"E:\pythonproject\simulation\HIRF\单波形分解.txt","w",encoding="utf-8")
    for row in range(len(wave_single)):
            file.write(f"{wave_single[row]}\n")
    file.close()

    evaluate_2wave=evaluate[0:5,:]
    evaluate_4wave=evaluate[5:10,:]
    evaluate_3wave=evaluate[10:15,:]
    combinations = [("evaluate_2wave", "evaluate_3wave"), ("evaluate_2wave", "evaluate_4wave"), ("evaluate_3wave", "evaluate_4wave")]
    wave_combined_table=[["wave_name_1","level","wave_name_2","level"]]
    for arr1_name, arr2_name in combinations:
        arr1 = locals()[arr1_name]  # 获取对应数组
        
        print(arr1)
        arr2 = locals()[arr2_name]
        
        print(arr2)
        print(f"\n{arr1_name} 和 {arr2_name} 的或运算结果:")
        
        # 对每一行进行或运算
        for i in range(5):  # 遍历每一行
            for j in range(5):
                row_result = np.logical_or(arr1[i], arr2[j])
                print(f"第一个波形{i+1}第二个波形{j+1}: {row_result}")
                if np.all(row_result):
                    print(f"当前波形可以组合为{arr1_name}的{i+1}电平和{arr2_name}的{j+1}电平")
                    wave_combined_table.append([arr1_name,i+1,arr2_name,j+1])
                else:
                    print(f"当前波形不可以组合为{arr1_name}的{i+1}电平和{arr2_name}的{j+1}电平")
    print(wave_combined_table)
    file = open(r"E:\pythonproject\simulation\HIRF\电压波形组合.txt","w")
    for row in range(len(wave_combined_table)):
        file.write(f"{wave_combined_table[row]}\n")
    file.close()

def wave_combind_current(peak,content,energy,rate_of_rise,rate_of_rise_1MHz):
    num_list=[["1波1电平",100,9.46e-3,5e-1,7e7],
              ["1波2电平",250,2.37e-2,3.13,1.75e8],
              ["1波3电平",600,5.68e-2,18,4.2e8],
              ["1波4电平",1500,0.142,112,1.05e9],
              ["1波5电平",3200,0.303,512,2.24e9],

              ["5A波1电平",150,1.96e-2,1.89,1.09e7],
              ["5A波2电平",400,5.22e-2,13.4,2.9e7],
              ["5A波3电平",1000,1.3e-1,84,7.24e7],
              ["5A波4电平",2000,0.261,336,1.45e8],
              ["5A波5电平",5000,0.652,2100,3.62e8],

              ["5B波1电平",150,1.02e-1,8.14,1.3e7],
              ["5B波2电平",400,2.72e-1,57.9,3.46e7],
              ["5B波3电平",1000,6.81e-1,362,8.66e7],
              ["5B波4电平",2000,1.36,1.45e3,1.73e8],
              ["5B波5电平",5000,3.41,9.04e3,4.33e8],


              ["6波1电平",5,0,0,6.65e8],
              ["6波2电平",12.5,0,0,1.66e9],
              ["6波3电平",30,0,0,3.99e9],
              ["6波4电平",75,0,0,9.98e9],
              ["6波5电平",160,0,0,2.13e10],


              ["3波1电平",20,0,0,1.33e8],
              ["3波2电平",50,0,0,3.33e8],
              ["3波3电平",120,0,0,7.98e8],
              ["3波4电平",300,0,0,2e9],
              ["3波5电平",640,0,0,4.26e9],

              ]
    
    # num_list_2wave=[
    #           ["2波1电平",50,1.18e-4,1.96e-3,2.73e8],
    #           ["2波2电平",125,2.94e-4,1.22e-2,6.83e8],
    #           ["2波3电平",300,7.07e-4,7.05e-2,1.64e9],
    #           ["2波4电平",750,1.77e-3,4.41e-1,4.1e9],
    #           ["2波5电平",1600,3.77e-3,2,8.74e9],
    # ]
    # num_list_4wave=[
    #           ["4波1电平",50,4.73e-3,1.25e-1,3.5e7],
    #           ["4波2电平",125,1.18e-2,7.81e-1,8.75e7],
    #           ["4波3电平",300,2.84e-2,4.5,2.1e8],
    #           ["4波4电平",750,7.09e-2,28.1,5.25e8],
    #           ["4波5电平",1600,1.51e-1,128,1.12e8],
    # ]
    # num_list_3wave=[
    #           ["3波1电平",100,0,0,6.65e8],
    #           ["3波2电平",250,0,0,1.66e9],
    #           ["3波3电平",600,0,0,3.99e9],
    #           ["3波4电平",1500,0,0,9.98e9],
    #           ["3波5电平",3200,0,0,2.13e10],
    # ]
    evaluate=np.zeros([15,4])
    for i in range(len(num_list)):
        if num_list[i][1]>peak:
            evaluate[i][0]=True

    for i in range(len(num_list)):
        if num_list[i][2]>content:
            evaluate[i][1]=True

    for i in range(len(num_list)):
        if num_list[i][3]>energy:
            evaluate[i][2]=True

    for i in range(len(num_list)-5):
        if num_list[i][4]>rate_of_rise:
            evaluate[i][3]=True

    for i in range(len(num_list)-10):
        if num_list[i+10][4]>rate_of_rise_1MHz:
            evaluate[i+10][3]=True
    print(evaluate)
    wave_single=[["wave_name_1&level"]]
    for row in range(len(evaluate)):
        if np.all(evaluate[row]):
            print(f"当前波形可由{num_list[row][0]}描述")
            wave_single.append({num_list[row][0]})

    print(wave_single)
    file = open(r"E:\pythonproject\simulation\HIRF\单波形分解.txt","w",encoding="utf-8")
    for row in range(len(wave_single)):
            file.write(f"{wave_single[row]}\n")
    file.close()

    evaluate_2wave=evaluate[0:5,:]
    evaluate_4wave=evaluate[5:10,:]
    evaluate_3wave=evaluate[10:15,:]
    combinations = [("evaluate_2wave", "evaluate_3wave"), ("evaluate_2wave", "evaluate_4wave"), ("evaluate_3wave", "evaluate_4wave")]
    wave_combined_table=[["wave_name_1","level","wave_name_2","level"]]
    for arr1_name, arr2_name in combinations:
        arr1 = locals()[arr1_name]  # 获取对应数组
        
        print(arr1)
        arr2 = locals()[arr2_name]
        
        print(arr2)
        print(f"\n{arr1_name} 和 {arr2_name} 的或运算结果:")
        
        # 对每一行进行或运算
        for i in range(5):  # 遍历每一行
            for j in range(5):
                row_result = np.logical_or(arr1[i], arr2[j])
                print(f"第一个波形{i+1}第二个波形{j+1}: {row_result}")
                if np.all(row_result):
                    print(f"当前波形可以组合为{arr1_name}的{i+1}电平和{arr2_name}的{j+1}电平")
                    wave_combined_table.append([arr1_name,i+1,arr2_name,j+1])
                else:
                    print(f"当前波形不可以组合为{arr1_name}的{i+1}电平和{arr2_name}的{j+1}电平")
    print(wave_combined_table)
    file = open(r"E:\pythonproject\simulation\HIRF\电压波形组合.txt","w")
    for row in range(len(wave_combined_table)):
        file.write(f"{wave_combined_table[row]}\n")
    file.close()

if  __name__=="__main__":
    # x=np.loadtxt(r"E:\pythonproject\HIRF\辐射敏感测试数据.txt")
    # x=pd.read_csv(r"E:\pythonproject\simulation\HIRF\波形分解资料\示例波形\tek1437.csv",skiprows=21)
    x = np.loadtxt(open(r"E:\pythonproject\simulation\HIRF\波形分解资料\示例波形\tek1437.csv","rb"),delimiter=",",skiprows=21)
    # print(x)
    time,amplitude=data_prepocessing(x,5000,200000)
    # print(time[1]-time[0],1/(time[1]-time[0]))
    # print("kkkkkk"*10)
    combined =zip(time,amplitude)
    file = open(r"E:\pythonproject\simulation\HIRF\lightning_raw.txt","w")
    for data in combined:
        file.write(f"{data[0]} {data[1]}\n")
    file.close()
    peak,content,energy,rate_of_rise,rate_of_rise_1MHz=lightning_abstract_separate(time,amplitude)
    print(f"当前峰值为{peak},电荷量为{content},能量为{energy},上升沿为{rate_of_rise},1MHz分量上升沿为{rate_of_rise_1MHz}")
    wave_combind_voltage(peak,content,energy,rate_of_rise,rate_of_rise_1MHz)
    
        

    
